﻿using System;
using System.Collections.Generic;
using Allegro.MathInterface;

namespace Allegro.Integrator
{
    public class SimpsonIntegrator : BaseIntegrator, IIterativeFunctionIntegrator
    {
        //const double EPS = 1.0e-6;
        const int JMAX = 20;

        /// <summary>
        /// Returns the integral of the function f from a to b. The parameters EPS can be set to the 
        /// desired fractional accuracy and JMAX so that 2 to the power JMAX-1 is the max. allowed
        /// number of stpes. Integration is performed by Simpson's rule
        /// </summary>
        /// <param name="f"></param>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <returns></returns>
        public double Integrate(Func<double, double> f, double a, double b)
        {
            double s, st, ost = 0.0, os = 0.0;
            
            for(int j = 1; j <= JMAX; j++)
            {
                st = Trapezoid(f, a, b, j);
                s = (4.0*st - ost)/3.0;
                if(j > 5)
                {
                    if(Math.Abs(s - os) < EPS * Math.Abs(os) || (s == 0.0 && os == 0.0))
                        return s;
                }
                os = s;
                ost = st;
            }
            throw new ApplicationException("Too many stpes in routine SimpsonIntegrator.Integrate");
        }
    }
}
